#clean up
rm(list=ls())

#preferelaes
options(scipen=8)

#load required libraries
library(spdep) 
library(maps)
library(foreign)
library(maptools)
library(spatstat)
library(mgcv)
library(ggmap)
library(geoR)
library(akima)
library(scatterplot3d)
library(sandwich) 
library(multiwayvcov) 
library(lmtest)
library(stringr)
library(xtable)
library(lattice)

####LOAD DATA###
#setwd("/Volumes/MONOGAN/registrationSpatial/LA/") ###SET TO YOUR WORKING DIRECTORY
la<-read.dta('laRegistration.dta')
la<-la[order(la$year,la$parish),]

###DESCRIPTIVES### 
summary(la)
boxplot(la$pctrep~la$year,xlab="Year",ylab="Proportion GOP Registration",ylim=c(0,.5))
cor(la$lnOil1,la$lnGas1)

#county list common to all maps
counties.to.map<-str_trim(paste('louisiana', str_replace(tolower(la$parish[la$year==1966]),pattern="[.]",replacement=""), sep=','))
counties.to.map[50]<-"louisiana,st martin:north" 
counties.to.map[65]<-"louisiana,st martin:south" 
counties.to.map<-sort(counties.to.map)
#counties.to.map

###OPEN MAP AND ENCODE SPATIAL CONNECTIONS FOR WEIGHT MATRIX###
#call la map
louisiana<-map("county", region=counties.to.map, fill=TRUE, plot=TRUE, resolution=0) 
IDs <- sapply(strsplit(louisiana$names, ","), function(x) x[2])
IDs[50:51]<-"st martin"
la.map.0<-map2SpatialPolygons(louisiana, IDs=IDs)

#create neighbor list from county map
la.nb<-poly2nb(la.map.0)
co.coords<-coordinates(la.map.0) #returns centroids
plot(la.nb, coords=co.coords)

###SPATIAL LAG FOR ALL WAVES USING KRONECKER PRODUCT###
wmat<-nb2mat(la.nb)
N<-64
Time<-49
ident<-diag(1,nrow=Time)
bigW<-(ident%x%wmat)
la$SpatLag<-bigW%*%log(la$pctrep)
la$LagPctblkreg<-bigW%*%log(la$pctblkreg)#lag of black registration

###PROBIT VERSION OF TABLE 6.2 WITH OIL, NO FIXED EFFECTS, AND SLX OF BLACK REGISTRATION###
prob.6.2.2<-lm(qnorm(pctrep)~qnorm(pctrep_t1)+ pctblkreg+ LagPctblkreg+pctblack +pctinmig+sla+pcturban+mhhia+lnOil1, data=la); summary(prob.6.2.2)

###LOG-LIKELIHOOD OF SPAITAL AND TIME LAG###
spatiotemporal <- function(par, X, y, W, Time, N) {
	n<-N*Time
	rho<-par[1]
	sigma2<-par[2]
	beta<-par[3:length(par)]
	ident<-diag(1,nrow=Time)
	nIdent<-diag(1,nrow=N)
	bigIdent<-ident%x%nIdent
	bigW<-(ident%x%W)
	first<-prod(1-rho*eigen(W)$values)^Time #Accurate and faster to use little W and exponentiate.
    loglikelihood <- log(first)-(n/2)*log(2*pi)-(n/2)*log(sigma2)-((t(y-rho*bigW%*%y-X%*%beta)%*%(y-rho*bigW%*%y-X%*%beta))/(2*sigma2))
    return(loglikelihood)
}

###SLX SPATIAL MODEL###
la.X<-as.matrix(cbind(1,subset(la,select=c(PROBpctrep_t1, pctblkreg, LagPctblkreg, pctblack,pctinmig,sla,pcturban,mhhia,lnOil1))))

la.spatiotemp <- optim(par=c(0,1,prob.6.2.2$coef),	# starting value for beta and sigma
   spatiotemporal,     									# the log-likelihood function
   method="BFGS",               						# optimization method
   hessian=TRUE,                						# return numerical Hessian
   control=list(fnscale=-1,trace=TRUE),    				# maximize function instead of minimize, trace progress
   y=qnorm(la$pctrep), X=la.X,W=wmat, N=64, Time=Time) 	#49 time points
   
#results   
la.spatiotemp$par
sqrt(diag(solve(-la.spatiotemp$hessian)))
la.space.time<-cbind(la.spatiotemp$par,sqrt(diag(solve(-la.spatiotemp$hessian))),la.spatiotemp$par/sqrt(diag(solve(-la.spatiotemp$hessian))),(1-pnorm(abs(la.spatiotemp$par/sqrt(diag(solve(-la.spatiotemp$hessian)))))))
la.space.time
colnames(la.space.time)<-c("Estimate","Std. Err.","z-ratio","p-value")
rownames(la.space.time)<-c("Spatial Lag","Error Variance","Intercept","Temporal Lag","Black Electoral Strength","Lag of Black Electoral Strength","Proportion Black","In-Migration","Southern Parish","Proportion Urban","Median Household Income","Logged Oil Production")
print(xtable(la.space.time,caption="Space-Time Model of Louisiana Parish-Level GOP Registration, 1966-2014",align="lrrrr",digits=4),type='html')

